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I. 


INTRODUCTION 


Our previous numerical simulations of the nonlinear evolution of 
the collisional Rayleigh-Taylor instability in the nighttime equatorial 
ionosphere (spread F) were confined to small (~ 3km) horizontal scale 
length initial perturbations and hence to furly developed spread F 
"bubbles" of approximately the same size in horizontal extent [ Scannapieco 
and Ossakow , 1976; Ossakow et al ., 1979], although spatially large 
vertically. However, observations by McClure et al [1977] also indi¬ 
cate ionospheric ion density "biteouts" of much larger horizontal extent 
(10 - > 200 km) and greater intensity (ion density depletions up to 
three orders of magnitude) than indicated by our small scale simulations. 
Therefore, we have extended our previous calculations and have performed 

a series of numerical simulations using much larger horizontal length scales 

\ 

(~ 75km) for our initial perturbations. These seed long wavelength 
perturbations, for example, could be due to neutral atmosphere gravity 
wave effects [ Rottger , 1976; Klostermeyer , 1978; Booker , 1979]. At 
the same time we have made very substantial improvements in the nu¬ 
merical techniques used to perform the simulations, including the 
utilization of the recently developed fully multidimensional flux- 
corrected transport (FCT) techniques of Zalesak [1979] . The results 
of our simulations indicate the following: 

1) large horizontal scale length initial perturbations evolve 
nonlinearly into large horizontal scale length equatorial 
spread F "bubbles;" 

2) these bubbles evolve on approximately the same time scale 
as do their smaller horizontal scale length counterparts; 


Note: Manuscript submitted December 4, 1979. 
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3) the plasma comprising these bubbles has its origin at much 
lower altitudes than that of the smaller horizontal scale 


length bubbles, resulting in plasma density depletions of 
very close to 100%. 

This last result is due to the fact that the polarization fields 
produced by ionospheric irregularities, aligned vertically, possess a 
fringe field component whose vertical extent is proportional to the 
horizontal extent of the irregularities producing the field. This is 
quite similar in origin to the fringe field produced at the edge of 
a parallel plate capacitor. Since the vertical extent of this fringe 
field determines the minimum altitude from which the rising bubble can 
draw plasma, it is not surprising that larger horizontal scale bubbles 
are more severely depleted. In Section II we give a brief review of the 
relevant theory, and of the basic equations used in our simulations. 
Section III contains the numerical simulations and a physical inter¬ 
pretation of the results is given. Section IV contains a summary, and 
in the appendix we describe briefly the numerical techniques used in our 
present computer code, emphasizing the differences between the present 
code and our previous one. 
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II. THEORY 


The geometry of the physical problem we are modeling is the same 
as in Ossakow et al . [1979]. All our simulations are carried out in 
a two dimensional (x,y) coordinate system. The constant magnetic field 
JJ is aligned along the z axis (pointing north). Gravity is directed 
along the negative y axis. Since our equations are two dimensional, 
n(y), V R (y), and V^(y) are taken to be representative values of the 
ambient electron density, recombination coefficient, and ion-neutral 
collision frequency (the result of integrating these quantities along 
magnetic field lines and dividing by a normalizing scale length). Mag¬ 
netic field lines are assumed to terminate at both ends in an electrically 
insulated medium (currents must close in the two dimensional plane, not 
in some distant E region). 

Following Ossakow et al . [1979], we describe the system with the 
two-fluid plasma continuity and momentum equations: 


St 
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(n v ) =* - v_n 
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where the subscript a denotes the species (i for ions, e for electrons), 
n is the species number density, v is velocity, V D is the recombination 
coefficient, E is the electric field, £ is the gravitational acceleration, 
q is the species charge, v is the species collision frequency with the 










neutral atmosphere, U is the neutral wind velocity, c is the speed 
—n 

of light, and m is the species mass. 

Note that, in contrast to Ossakow et al [1979], we have dropped 
the term + V^n^ from (1). This is the equivalent of dropping the 
assumption of the existence of an ionization source given by that 
term. This ionization source was such that the ambient ionization 
profile n (y) was an equilibrium profile (dn_^/dt = 0). Our present 
model therefore has instead 


dn 

0.0 

- ■ = - v„n 

dt R ao 


(3) 


Hence, when normalized results n^(x,y)/ n ao (y) are later presented, 
it should be understood that both the numerator and denominator are 
time dependent. 

Figure 1 shows the recombination rate V R and ion-neutral collision 

frequency used in our simulations. It shall be seen presently 

that V gn need not be specified as long as it is much smaller than 

the electron gyro frequency ft g . For details on the form of and 

V R as depicted in Figure 1, see Ossakow et al .[l979]. If we neglect 

the inertial terms (the left-hand side) of (2) by assuming the inertial 

time scales are much larger than either the gyro periods or the mean 

time between collisions, then the equation can be inverted to give an 

algebraic expressions for v^. In two-dimensional (x,y) geometry with 

15 along the z axis, the solution is for our problem, with m g « , 

V. Al, « 1, V /ft =0 (where ft = eB/m c), and U =0. 
in i en e a <X ~n 
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We now make the electrostatic approximation, 


E = V $ 


( 6 ) 


where V * x(a/dx) + y(3/Sy)»and the quasi-neutrality approximation 

n ~ n s n • We then have 
e i 


V x . j = 0 (7) 

pen - vj (8) 


Substituting (4) and (5) into (8) and evaluating (7), we have for 
the electrostatic potential: 


(V in n V ) = 


” e 8 3y (V in n) ' - 8 dx 


(9) 


As in Ossakow et al . [1979] we set $ = where V x ^ Q = 

A 2 

- (m^g/e) y. Since 7^ 4> o = 0, our final potential equation becomes 
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The effect of <f> Q is merely to superimpose a bulk westward plasma 
velocity on the electron velocity field determined from $ , 

without affecting the morphology of the developing structures. 

Hence, we ignore this motion. 

Our assumption of quasi-neutrality has made one of our two 
continuity equations (1) redundant. We therefore choose the electron 
equation for its simplicity: 


dn d / nc 

3t dx \B 


^l\ 

ay / 




^i\ 

dx ) 


( 11 ) 


III. NUMERICAL SIMULATION RESULTS AND DISCUSSION 

Equations (10) and (11), together with appropriate boundary 
conditions, constitute the nonlinear system of equations we shall solve 
numerically. Note that in contrast to Ossakow et al . [1979], we 
have chosen not to put the equations into a normalized form by dividing 
through by n Q (y). The numerical techniques used to solve these equa¬ 
tions are discussed in the appendix. 

The numerical calculations to be presented were performed on 
a two-dimensional cartesian (x,y) mesh using 42 points in the x (east- 
west) direction, and 142 points in the y (vertical)direction. The 
(uniform) grid spacing was 2km in the y direction for all calcu¬ 
lations. The grid spacing in the x direction was 200m in the "small" 
horizontal scale length cases and 5km in the "large" cases. The 
bottom of the grid corresponds to 252km altitude and the top of the 
grid to 534km altitude in all simulations. Periodic boundary con- 


6 




ditions were imposed on both n and 4>^ in the x-direction. In the y 
direction transmittive boundary conditions were imposed on n(dn/dy ■ 0) 
and Neumann (d<^/dy = 0) boundary conditions were imposed on 

Three kinds of plots will be presented: (1) contours of con¬ 
stant n(x,y,t); (2) contours of constant n(x,y,t)/n (y,t); and (3) con- 

o 

tours of constant electrostatic potential 4>^. Superimposed on each 

contour plot is a dashed line depicting n (y.t) for reference purposes. 

o 

Our n o (y,o) profile is such that the F2 peak is located at 434km alti¬ 
tude, and the minimum electron density scale length L=n (dn/dy) 3 

“ o o 

is 10km. The simulations are all identical except for the east-west 
grid spacing Ax and the form of the initial perturbation. Two kinds 
of initial perturbations were used: 

= 1 - e' 3 cos , 0 £|x|£ 8Ax 

■ 1 } [-(Sr) - 1 ] • 

8 Ax £|x|^ 16Ax 

= 1 , | x | > 16 Ax 

( 12 ) 

• 1 - e " 3 c ° s (ifc) <13 > 


n(x,y,0) 
n (y,0) 


Perturbation A: 


n(x,y,0) 
n (y,0) 


n(x,y,0 ) 
n (y,0) 


Perturbation B: 


n(x,y,0) 
n (y,0) 
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Perturbation A is exactly the form used in Ossakow et al. [1979], 
and perturbation B is a pure sine wave of wavelength AO Ax (our system 
length in the x direction). Both represent maximum initial pertur¬ 
bation amplitudes of approximately 5%. Four simulations have been run: 
(i) lS-Perturbation A with Ax = 200m; (ii) lL-Perturbation A with 
Ax = 5km; (iii) 2S-Perturbation B with Ax = 200 m; and (iv) 2L-Per- 
turbation B with Ax = 5km. Calculation IS above is identical to 
ESF III of Ossakow et al . [1979]. The "large" versus "small" 
comparison obviously involves comparing calculation IS to calculation 
1L, and calculation 2S to calculation 2L. One notes that for the 
minimum L 10km in our simulation, kL > 1 for the IS and 2S cases 
and kL < 1 for the 1L and 2L cases. 

Figure 2 shows isodensity contours of calculation IS at times 
300, 700, 1000, and 1200 seconds after initialization. Figure 3 shows 
the same contours at the same times but for calculation 1L. The 
presence of much lower density fluid in the bubble in calculation 1L 
is obvious. Also obvious is a basic difference in the bubble mor¬ 
phology at late times. At 1200 seconds, iS has pinched off into two 
bubbles, with the more intense one below the initial central bubble. 

In addition, another bubble has formed in the sides of the calcu¬ 
lation . These structures are more obvious in the plot of n(x,y)/n Q (y) 
at 1200 seconds for IS shown in Fig. Aa. The maximum depletion levels 
are 70% in the top central bubble ,97% in the lower central bubble, 
and 95% in the side bubble. Note that here and in all subsequent 
plots of (n,x,z)/n Q (z) the contour plotting is such that the first 


8 




(outer) depletion contour n/n^ is 0.5 and each succeeding inner 

contour is 0.5 times the previous one. For example, the lower bubble 

in Fig. 4a has five contours. The outermost would have n/n = 0.50 

o 

(50% depletion), the next inner one n/n Q = 0.25 (75% depletion), the 
second inner one n/n Q = 0.125 (87.5% depletion) and the innermost 
n/n^ = 0.03 (97% depletion). For the enhancement contours (dashed 
lines) the first outer contour is 2.0 and the succeeding inner ones 
are 2.0 times the previous ones. In obtaining percentage enhance¬ 
ments and depletions one then subtracts 1.0. 

Calculation 1L, on the other hand, shows a single plume of 
depleted ionization at 1200 seconds, with no secondary central 
bubble and no side bubble. There also is no indication of a widen¬ 
ing of the top of the bubble, as there is in IS. In Fig. 4b we 
show a plot of n(x,y)/n Q (y) for 1L at 1200 sec. The level of 
depletion is greater than 99.9% for the entire 10km by 70km oval 
"hole" located inside the tenth solid contour of Fig. 4b and 
represents at least a three order of magnitude decrease (biteout) 
in plasma density. 

We now repeat the above comparison, but this time for per¬ 
turbation B (calculations 2S and 2L). Figure 5 shows isodensity 
contours of calculation 2S at times 300, 700, 1000, and 1091 
seconds after initialization; while Figure 6 shows similar plots 
of calculation 2L at times 700, 1000, 1200 and 1364 sec. Comparison 
again shows the presence of much lower density plasma in the bubble 
in the 2L calculation. Morphological differences are also present, 
the most notable being the widening of the top of the bubble in 2S 
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which is not present in 2L. Figure 7 shows a comparison of the n/n 

o 

profiles at late time. Again maximum depletions in the 2S case are 
about 97%, while a large portion of the 2L plume is 99.9% depleted or 
greater. 

We can also compare the effect of the form of the perturbation 
by comparing IS to 2S and 1L to 2L. The latter comparison shows striking 
similarity, whereas the former shows some differences, the most notable 
being the lack of central bubble "pinching" and the lack of lateral 
bubbles in case 2S. We conclude that the morphology of the late-time 
bubbles is dependent, at least somewhat, on the form of the initial per¬ 
turbation. 

Bubble rise velocities are of some interest, and we give below the 
average bubble rise velocity for each case, computed from the last two 
frames of Figs. 2, 3, 5, and 6: 


IS 

210 m/sec 

1L 

230 m/sec 

2S 

420 m/sec 

2L 

280 m/sec 


The rise velocity of an individual bubble is dependent upon the rela¬ 
tive depletion level of the bubble, its geometry, and upon interactions 
with other plasma structure nearby. These first two effects are treated 
in Ossakow and Chaturvedi [1978], and the present results above are con¬ 
sistent with the results therein. For instance, the relatively high 
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rise velocity associated with 2S is seen to be due to the fact that 
the bubble is more severely depleted than that in IS. Further, the 
roughly equal rise velocities of the IS, 1L, and 2L bubbles, in spite 
of the fact that the 1L and 2L bubbles are much more severely depleted 
than that in IS, is explained by noting that IS actually approximates 
the geometry of a "sheet" of depleted plasma, whereas the 1L and 2L 
bubbles more closely resemble a cylindrical geometry (see Ossakow 
and Chaturvedi [1978]). 

In an attempt to understand the reasons for the differences 
in the nonlinear evolution of small and large horizontal scale per¬ 
turbations, we look at the potential equation: 


V 2 


V(V in n ) 

Vin n 


V <p 


= 

1 C V 


in 


1 dn 
n dx 


At early times we expect V to be small with respect to Bg/cv^, 
so we ignore the second term on the left hand side, giving a Poisson 
equation for <j>^: 

2 A _ Bg 1 dn 

V ^1 cV. n dx 

in 

We can now interpret the right hand side as simply the local charge 

density p (such that V.JE = p) . Since we have initialized all of these 
1 

calculations with — — independent of y (see (12)), what we are deal¬ 
ing with is a distribution of charge density that has the form of 
diffuse "plates" aligned in the vertical direction. Noting that the 
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term v. decreases with altitude, we find that these diffuse "plates" 
in 

have an equally diffuse "edge" in the y direction. Taking the analogy 
to its conclusion, we model our initial conditions, or any vertically 
aligned structure for that matter, as an array of plates of charge 
(non-conducting capacitor plates) with an edge somewhere in or above 
our grid. Tn Figure 8 we show schematically the electric potential 
field surrounding the edge of parallel plates of charge for two 
different separation distances. Since there is only one scale length 
in the configuration (the plate separation distance d), then all other 
scale lengths must be proportional to this distance. In particular, 
the characteristic distance parallel to the plates over which the 
electric field outside the plates (the fringe field) falls off must be 
proportional to d. Since in our problem the contours of electrostatic 
potential are in fact streamlines (see (11)), this distance will 
determine the maximum depth in the fluid from which the electrostatic 
field will draw fluid into the bubble. Since the plasma density is 
lower at greater depths, this distance will determine the minimum 
plasma density inside the bubble. To test these ideas, we examine 
the actual early time electrostatic potential fields from the cal¬ 
culations we have nresented. Fig. 9 shows contours of n and 
for the IS initial conditions, and the same plots for 1L. A similar 
comparison for cases 2S and 2L is shown in Fig. 10. All contour 
plots of are scaled in such a way as to evenly space exactly 12 
contour lines between the maximum and minimum value of , to 
normalize the plots so they can be compared without bias. The 
increased vertical extent of the contours of (streamlines) for 
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cases 1L and 2L are evident. 

Of course, the initial profile generating these aligned plates 
of charge lasts only a short time. In the linear phase of growth, 
the perturbation grows in the region of linear instability (the F region 
bottomside), and damps elsewhere. Our "plates" therefore very quickly 
become horizontally spaced regions of charge with a limited vertical 
extent, confined to the region of steepest vertical gradient on the F 
region bottomside. Nonetheless, the scaling arguments advanced above 
still hold: the vertical extent of the polarization electric field 
scales as the horizontal extent of the structure causing the field. 

This is easily seen in Figure 11 and 12 where comparison is made of the 
<j>l contours for IS vs 1L and 2S vs 2L respectively, at a time of 700 
sec. Cases IS and 2S are seen to be mixing fluid over a fairly narrow 
altitude range, while 1L and 2L have each formed a large convective 
cell more than 150km in vertical extent, drawing plasma fluid into the 
bubble from deep in the ionosphere. It is not surprising, then, that 
the larger horizontal scale bubbles are more severely depleted at 
late times. 

IV. SUMMARY 

On the basis of our numerical simulations, and of a qualitative 
scale analysis of the driving electrostatic potential equation, we 
conclude that the severe "biteouts" of three orders of magnitude and 
bubble rise velocities of 150 m/sec reported by McClure et al [1977] 
are completely consistent with large (~75-200km) horizontal bubble 
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size scales. In our simulations, the severe biteouts associated 
with large horizontal scale lengths are due to the fact that the 
plasma comprising these bubbles has its origin at much lower alti¬ 
tudes then in the small horizontal length scale cases. Again, these 
results are consistent with those of McClure et al. [1977], who 
base their conclusions on ion mass spectrometer measurements of the 
H + - 0 + - N + balance inside the bubbles, which they find to be 
'^characteristic of undisturbed plasma found at lower altitudes." 

The variation in the vertical velocities of various bubbles noted 
by McClure et al. Ll977] is probably due to interactions between 
bubbles. Note, for example, that in Fig. 2, the secondary bubble 
is rising at a much slower rate than is the central bubble. Bubble 
interaction will be the subject of forthcoming theoretical and 
numerical studies. 
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APPENDIX: Numerical Solution of the Equations 

Of the two partial differential equations we must solve, 

(10) is elliptic and linear and (11) is hyperbolic and nonlinear. 

Both equations represent numerical challenges, and we could easily 
devote the bulk of this paper to the numerical techniques used for 
their solution. However, as we stated in the introduction, we shall 
confine ourselves to the improvements made in these techniques since 
the calculations of Ossakow et al . [1979]. We begin with (11). 

Equation (11) is solved in finite difference form 
using a fully multidimensional second order in time, fourth order in 
space, leapfrog-trapezoidal flux-corrected transport (FCT) scheme. 

Both the higher order leapfrog-trapezoidal scheme itself, as well 
as the fully multidimensional algorithm utilized in the critical 
flux-limiting stage of FCT, are recent developments and are described 
by Zalesak [1979]. These developments represent significant extensions 
of the theory of FCT, a numerical technique originated by Boris and 
Book [1973] to handle equations of the form (11) where steep gradients 
are expected for form. By contrast, the calculations in Ossakow et al. 
[1979] used a first order in time, second order in space FCT algorithm 
which was only one-dimensional (since fully multidimensional FCT 
algorithms did not exist at the time), and hence used time-splitting 
(sequential x and y operators) to solve the two-dimensional equation 

(11) . It is known that time-splitting can introduce numerical 
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problems into an incompressible flow calculation (see Zalesak , 1979), 
although our previous equatorial spread F calculations did not 
exhibit any of the symptoms of these difficulties. 

Equation (10) is the elliptic equation whose solution 
gives us the electrostatic potential $ . The right hand side is 
known and the left hand side represents a Hermitian operator 
operating on 4>^, giving only real eigenvalues and apparently easing 
the difficulty of solution. However, the extremely large range for 
the values of the quantity n makes for an equally large span of 
operator eigenvalues, and solution of the equation in this form 
using iterative techniques has not been successful. We have found 
one (and only one) method of direct solution, the stabilized error 
vector propagation (SEVP) method of Madala [1978], but the execution 
speeds for SEVP are not as favorable as for the method we now describe. 

We start by expanding the operator and dividing through 
by v. n, as was done by Ossakow et al . [1979], giving 


V(v, n) 
in 


V $ = - 

a 1 c V. n dx 
in 


The equation is now in a form suitable for solution by the Chebychev- 
iterative relaxation technique of McDonald [1977]. However, great 

care must be given to the differencing of the term ~ - V(v^ n n) and 

1 ^ 1 ^ 

— -r— ,and this is the point we wish to address. We work with the 
n 3x 

term ^ in one spatial dimension, since this example will make 

the point. A straightforward second-order difference form for this 
term is 
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(A-2) 



where the subscript i refers to grid point location in the x direction 
and Ax is the (uniform) spacing between grid points. This is the 
form used in Ossakow et al . [1979J. We shall show that this difference 
form produces solutions <}>^ with potentially undesirable properties, 

and causes undue numerical hardship in finding those solutions. 

1 dn B 

Let us rewrite — as (In n). A physical inter- 

n ox ox 

pretation of the term is now much easier: the term tells us how 
rapidly the logarithm of n is varying with respect to x. Suppose, 
for argument's sake, that the smallest and largest values of n in 
the problem are 10^ and 10^ respectively. On a grid of size Ax', the 
largest value representable for that term would occur when a fluid 

element of density 10^ and one of density 10^ occupy adjacent grid 

£ 

points. The value of — (In n) evaluated midway between these two 
ox 

grid points would be ^ ( In 10^ - In 10 1 ) = 9.2 Ax Evaluation 
of (A-2) for n^ ^’"i’ an< * n i+l ^ av ^ n 8 values of 10 , 10^, and 10^ 
respectively gives a value for (^ -^) of 5 x 10^ Ax \ far in 
excess of the maximum value for this term given by the above argu¬ 
ment. Logarithmic interpretation of this term would state that n 
4 

varied by more than 10 orders of magnitude over a single grid 
spacing, a ridiculous statement in light of the fact that there 
are only four orders of magnitude of n in the problem. 
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of the term 


The potentially large values of these terms, in particular 


- V (v, n) in (A-l), not only cause extremely slow 

n j.' in 
in 

convergence of the iterative solution, but can also put spurious 


oscillations into the exact finite difference solution due to 
cell Reynolds number effects [ Roache , 1976]. As shown by Roache 


[1976] these oscillations can occur any time the value of the 

term —-— V (V. n) in (A-l) exceeds a critical value of 2 Ax 
V. n jl in 
in 

It is clear, then, that (A-2) represents an undesirable difference 


form for these logarithmic terms. Better approximations are 


,1 dn. 1 

n dz . 2Ax n 


'i+1 - ln n i-l } 


(A-3) 


and 


(i 



1 _ 

Ax 


"i+l ~ n i-l 
n i+l +n i-l 


(A-4) 


Equation (A-3) is probably the most accurate, but evaluation of 
the logarithms at every time step is computationally expensive, and 
for problems like the one at hand where n varies by many orders of 
magnitude across the grid, there is still no guarantee that the 
critical cell Reynolds number will not be exceeded. For these 
reasons we use (A-4) for the present calculations. 
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Fig. 2 - Sequence of four plots showing iso-electron density 
contours of calculation IS at 300, 700, 1000, and 1200 sec. 
Superimposed on each plot is a dashed line depicting n (z,t). 
Electron densities are given in cm . The observer is looking 
southward. 
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